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We show that the dynamical melting of a Mott insulator in a three-dimensional lattice leads to 
condensation at nonzero momenta, a phenomenon that can be used to generate strongly interacting 
atom lasers in optical lattices. For infinite onsite repulsion, the case considered here, the momenta 
at which bosons condense is determined analytically and found to have a simple dependence on the 
hopping amplitudes. The occupation of the condensates is shown to scale linearly with the total 
number of atoms in the initial Mott insulator. Our results are obtained using a Gutzwiller-type 
mean- field approach, gauged against exact diagonalization solutions of small systems. 
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The invention of the first optical lasers more than half a 
century ago marked the beginning of the ultimate control 
over light waves in a way that would later revolutionize 
the world. Now, with the realization of Bose-Einstein 
condensation (BEC) in atomic gases we are entering 
an era in which a similar degree of control over matter 
waves is being achieved. In recent years, much effort has 
been devoted to developing techniques for converting the 
trapped atoms of a BEC into freely propagating coherent 
matter waves - the so called atom lasers Q, Q . The real- 
ization of atom lasers [3[ has opened promising avenues 
of theoretical and experimental research. These coherent 
matter- wave beams are expected to become constituents 
in future scientific and technological devices 

So far, most of the atom-laser experiments have dealt 
with weakly-interacting gases. The possibility of gain- 
ing further control by enhancing interactions by means 
of optical lattices has remained open. Here we address 
this issue and demonstrate that the free expansion of 
bosons from a Mott-insulating state into an empty lat- 
tice leads to the emergence of coherent matter waves at 
nonzero momenta, which in turn can be fully controlled 
by the hopping amplitudes in the lattice. In one dimen- 
sion (ID), a related phenomenon was reported in Ref. [H|, 
where the melting of a Mott insulator of hard-core bosons 
(HCBs) was shown to produce a quasi-coherent matter 
wave characterized by power-law decaying correlations. 
Interestingly, the power law was found to be the same 
as for the ground-state solution, i.e., in this transient 
regime, the system cannot be characterized by a finite 
temperature. Analytical work in the XY chain repro- 
duced this behavior [6[, and further studies in ID showed 
that soft-core bosons Q and two-component fermionic 
systems [i[ exhibit similar phenomena. 

Following the promising results obtained for ID sys- 
tems, in this Letter we examine the dynamics of bosons 
in higher dimensions. For concreteness, we will focus on 
the hard-core limit ^ . Studying the expansion in higher 
dimensions is of much interest not only because most 
experiments are performed in these regimes but also be- 
cause true condensation, and hence a true atom laser, can 



only be realized in dimensions higher than one. We ex- 
amine whether condensation (off-diagonal long-range or- 
der) can develop dynamically during the expansion of an 
initial state that exhibits only short-range correlations, 
namely, a Mott-insulating state. 

In this study, we consider HCBs on a three dimensional 
(3D) lattice, with N = L x x L\ sites in the presence of 
a general onsite potential. Here, L x (L±) denotes the 
number of sites in the x (y and z) direction(s). The 
Hamiltonian can be written as: 

H = - ^2 Uj + a]a^j - ^ \nhi , (1) 

(ij) 1 

where {ij) denotes nearest neighbors, di (a\) destroys 
(creates) a HCB on site z, hi = d\di is the local density 
operator, \±i is the local chemical potential, and the t^-'s 
are the hopping amplitudes. The HCB creation and an- 
nihilation operators satisfy the constraints a\ 2 = af = 0, 
which preclude multiple site occupancies. Equation (pQ) 
may be viewed as a spin- 1/2 XY model with a site- 
dependent magnetic field applied in the z direction. 

A typical initial state in our system is an almost perfect 
Mott insulator, created by placing HCBs in a very strong 
harmonic trap with curvature along the x direction, i.e., 
with fii = fio — Vxf, where V is the curvature of the 
trap, X{ is the distance of the i-th site from the center 
of the trap along the x direction, and jao is the chemical 
potential at the center of the trap. The resulting initial 
state is a 'slab' of width W x around the center of the 
trap in the x direction with (hi) « 1, which extends 
over the entire lattice in the transverse directions. In 
experiments, such a state could be produced by a box- 
like trap in the transverse directions and very strong on- 
site interactions (U/t ^> 1). For simplicity, we consider 
periodic boundaries in the transverse directions. This 
choice makes our system homogeneous in those directions 
as a result of which, the actual value of L± (L± > 2) 
plays no role in the mean- field analysis. 

At time r = 0, the harmonic trap is turned off and 
the bosons start expanding freely throughout the lat- 
tice. Their particular initial state causes them to move 
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in the x direction leaving the system translationally in- 
variant in the transverse ones. We then study various 
properties of the system measured in the course of the 
evolution. Among those are density profiles, momentum 
distribution functions, and the condensate fraction along 
with the lowest natural orbital (LNO). The latter two are 
the largest eigenvalue of the one-particle density matrix 



e density i 
ively Ho(|. 



Pij = (a]dj) and its eigenvector, respect 

Partial control over the dynamics of the system is 
gained by considering two different hopping parameters: 
The first is the hopping amplitude in the x direction, de- 
noted by t x . Without loss of generality, this amplitude 
is fixed at t x = 1, setting the energy scale. The second 
value corresponds to the hopping amplitude in the trans- 
verse directions t y = t z . For reasons that will become 
apparent later, we shall use the dimensionless quantity 
rj = (t y + t z )/t x as the control parameter. 

In contrast to the special ID case, where the model 
has an exact solution |5], to explore the dynamics of 
HCBs in higher dimensions one must resort to approxi- 
mate schemes. This is because numerically-exact meth- 
ods such as exact diagonalization, scale exponentially 
with the system size and hence allow insight into the 
dynamics of only small systems. Here, we probe the dy- 
namics of HCBs primarily by means of a Gutzwiller-type 
mean- field approximation, adjusted to handle time evolu- 
tion via the time-dependent variational principle method, 
originally developed for soft-core bosons 1111. This ap- 
proach becomes exact in the limit of infinite dimensions, 
but has proved to provide qualitatively correct phase di- 
agrams for the ground state of two-dimensional (2D) and 
3D hard-core boson systems [l2| . 

Within this method, we employ the following simple 
product ansatz as the wave-function of the system: 
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where we have utilized the usual parametrization of spin- 
1/2 particles. The polar angles Oi G [0,7r], and the az- 
imuthal angles, Xi an d <j>i G [0,27r), are time-dependent. 
By stationarizing (ip\id / dr — H\ip) (here, h = 1), we ob- 
tain the following set of first-order coupled equations: 



0i = 2 Im [Si] 



= ^ - 2 Re [Si] cot Oi 



(3) 



where Si = \ £V Uj sin^e^' 0i) ^ n(j) , and S^^m is 
unity if i and j are neighbors and zero otherwise [13| . 

Equations ([3]) are then solved simultaneously for all 
sites using a multi-dimensional forth-order Runge-Kutta 
method with an 0(Ar 5 ) error, where At is the discrete 
time step (taken to be Ar = 0.001). Self consistency of 
the method is verified by monitoring the conservation of 
quantities such as total number of particles, total energy 
and total momentum. For the chosen time step, these 
quantities were observed to remain constant (up to 8 sig- 
nificant digits) throughout the entire simulation. 



Since our description of the model is approximate in 
nature, in order to gain an understanding of the accuracy 
of our results, we compare them against exact solutions. 
The latter are available in ID, where the mean- field ap- 
proach is not expected to be accurate [14], and for small 
2D systems where exact diagonalization calculations can 
be used. 

In Fig. [TJ we present an analysis of several small 2D 
systems. The figure shows the time evolution of the nor- 
malized difference for (nk) between the mean- field results 
and the exact ones, defined as 
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for lattices with up to 6 x 4 sites and up to six bosons, 
and two different values of r] [14[. 

As expected, the errors increase during the expansion 
and attain significant values in all cases (~ 0.8 in the 
worst case and « 0.15 in the best case) even before they 
start bouncing off the edges of the lattice (which occurs 
at r ~ 1). However, several conclusions may be drawn 
from comparison of the errors in the various cases. First, 
larger values of r] yield smaller errors [Fig. Ufa) vsQJb)]. 
This is not surprising considering that as rj increases, the 
system moves from ID (77 = 0) to 'full' 2D (77 = 1), and 
that mean- field theories are expected to become more ac- 
curate as the dimensionality of the system increases [l2[ . 
Figure [1] also shows that two other factors contribute to 
the reduction of the error, (i) adding particles to the sys- 
tem, and (ii) increasing the lattice size in the transverse 
direction 1141. 
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FIG. 1: (Color online) Momentum-distribution error 6mom 
[Eq.Q] as a function of time r for different small 2D systems 
of size L x x L± with Nt, bosons. The hopping ratios are: (a) 
77 = 0.2 and (b) 77 = 0.8 Q. 

The fact that in 2D, the mean- field results are qualita- 
tively similar to the ones obtained with exact diagonal- 
ization [Til , and the quantitative decrease of the relative 
errors as rj becomes large or as the number of particles 
increases, suggest that in 3D one can not only gain a 
qualitative understanding of the dynamics of the bosons 
during the melting of a Mott insulator, but also an actual 
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quantitative picture of it. This is both because mean-field 
is expected to be more accurate in 3D than in 2D, and 
because the systems we can simulate contain far more 
particles on much larger lattices. 

Figure [2] depicts the main result of this Letter. It shows 
the density profiles and momentum distribution functions 
of the evolution of A/5 = 40 x L\ bosons on a 700 x L^-site 
lattice, during the melting of a three-dimensional Mott 
insulator with 77 = 0.8. Several important features are 
apparent in these plots. In the density profiles (left pan- 
els), one can see that during the expansion, after the 
insulator has completely melted, two solitonic 'lumps' 
of bosons form and subsequently move in opposite di- 
rections with constant velocity. Even though the lumps 
slowly broaden as they travel, they retain their shape 
[compare Figs. [2je) and[2jg)]. Cross-sectional plots of 
Fig. [2] are presented in the supplementary material. 
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FIG. 2: (Color online) Evolution of: density profiles (left pan- 
els) and momentum distribution function, as a fraction of L\ 
with the z direction omitted (right panels), during the free 
expansion of 40 x L\ HCBs in a 3D lattice with L x — 700 
and 77 = 0.8. The times shown here are, from top to bottom, 
T = 0, r = 84, r = 151 and r 187. 



where |k) = aJjO) and (ak) creates (destroys) a HCB 
with momentum k. Note that even though the disper- 
sion relation is the same as for noninteracting particles, 
HCBs do interact due to the constraints on site occu- 
pancy. (This is reflected in Fig. [2j where momenta are 
redistributed during the expansion.) In our setup, the 
initial total energy of the system is e « 0, as can be 
immediately verified by considering the energy of the in- 
sulating ground state. Since the expansion takes place in 
the x direction, it is expected that if condensation occurs 
it would take place at momenta (=bg x ,0,0), for some q x 
which is to be determined. Further assuming that the 
majority of the bosons occupy these two modes, together 
with conservation of energy, leads [from Eq. (j5j)] to a sim- 
ple relation between the momentum of the condensate 
(q x ) and the parameters of the model: 

cosq x = -rj. (6) 

This remarkable behavior of the HCB gas therefore im- 
plies that the locations of the two prominent momentum 
modes are fully controlled by the hopping parameters. 

The mean-field results agree in full with this analyti- 
cal argument. This is shown in Fig. [3fa) where we plot 
the locations of the momentum peaks for different values 
of 77 as predicted by Eq. (j6j) (solid line) and the result 
of the mean-field simulation (points). The velocity at 
which the condensates travel in the lattice is then given 
by v g = Vek = 2t x s'mq x x. Interestingly, if the hopping 
parameters are such that the right-hand-side of Eq. (|6j) is 
greater than unity, the bosonic gas does not melt within 
the mean-field theory. This is an artifact of the mean- 
field approach. Exact-diagonalization calculations on 2D 
systems show that the Mott insulator does slowly melt for 
77 > 1. Hence, we expect that in large finite-dimensional 
systems with 77 > 1 the Mott insulator will melt but no 
condensation will occur. Since the 2D systems we can 
study with exact diagonalization are too small to ver- 



The above behavior of the density profiles is accompa- 
nied by an even more remarkable behavior of the momen- 
tum distribution function (right panels in Fig. [2j). While 
(h\g) of the initial state is almost completely flat (as the 
Mott insulator is localized in real space), as the system 
expands, interactions between bosons redistribute the en- 
ergy. This results in the emergence of two condensates at 
opposite momenta. Given the dispersion relation in the 
lattice, these momenta determine the velocity of the two 
lumps seen in the density profiles [14 1. 

As it turns out, the momenta at which condensation 
takes place can be evaluated analytically without resort- 
ing to mean-field theory. To calculate them, we first ob- 
tain the dispersion relation of the system in the absence 
of the trapping potential, given by 



ek 



(k\H\k) 



2t x cos k x 



2ty cos ky 



2t z cos k z (5) 



1 — 

♦ peak 

— cos -1 (-77) 


(a) J 






1- 



-80 ^ 



-40 



AO 







1 



50 100 
W x (width of initial state) 



0.5 

r] (hopping ratio) 

FIG. 3: (a) Dependence of the momentum peak locations on 
the hopping parameters. The data points are taken from the 
various simulations whereas the line is cos -1 (—77). (b) Scaling 
of the three-dimensional system: Total number of condensed 



bosons (in units of L ± ) as ; 
Mott insulating state for 77 
(x) and 77 = 0.8 (□). 



function of width of the initial 

= 0.2 (O), V = 0.4 (A), V = 0.6 



ify this hypothesis, this matter will have to be tested in 
experiments with ultracold bosons in optical lattices. 

So far, we have been stating without proof that during 
the expansion of the initial Mott insulator, the emergence 
of peaks at nonzero momenta signals condensation. This 
was certainly not the case in the ID systems studied in 
Ref. [5] , even though similar peaks were observed at mo- 
menta q x = ±7r/2. In ID, the occupation of the two 
LNOs [which correspond to the two peaks in n(k)] was 
found to be proportional to y/Nb, meaning that, in the 
thermodynamic limit, each condensate fraction (the ra- 
tio of the number of particles in each LNO to the total 
number of particles in the system) is zero True con- 
densation only takes place if the LNO occupation is in 
proportion to the total number of particles in the system 
[15j, hence one must verify that it scales with N^. 

In Fig. [3jb), we study the scaling of the total occupa- 
tion of the LNOs in our systems for different values of 
77 ranging from 0.2 to 0.8. As the figure indicates, the 
number of condensed bosons indeed grows linearly with 
the width W X1 i.e., with the total number of bosons in 
the initial Mott insulating state. As the value of rj in- 
creases and one moves from a quasi- ID system to a more 
3D system, the slope, that is, the fraction of bosons that 
condense, grows as well. While it is well known that 
true condensation can occur in 2D systems at zero tem- 
perature and in 3D or higher dimensional systems below 
some critical temperature, the surprising result here is 
that condensation is found out of equilibrium, in a sys- 
tem where the energy is conserved, and for an insulating 
initial state that exhibits no off-diagonal correlations. 

The dynamics of the systems analyzed above describes 
bosons hopping on a 3D lattice with periodic boundary 
conditions in the transverse directions. These were cho- 
sen in order to eliminate edge effects and to keep the 
system homogeneous in those directions. This geome- 
try is, unfortunately, not realizable experimentally. In 
experiments, the boundaries are naturally open (a box 
trap or other trapping potentials). In the supplementary 
material, we show that our conclusions are not altered 
when open boundaries (a box trap) are considered in 
the transverse directions provided that L±_ is sufficiently 
large. This is to be expected since L±_ —> 00 and peri- 
odic boundary conditions in the _L directions are equiv- 
alent. In pra ctice, the results are found to be similar for 
L± > 40 PJ, which are lattice sizes that are accessible 
experimentally. 

In summary, we have shown that the melting of a 
Mott insulator with (hi) ~ 1 in a three dimensional lat- 
tice leads to the dynamical emergence of condensates at 
nonzero momenta. We have determined the momenta at 
which condensation takes places and showed that these 
values have a simple dependence on the hopping ampli- 
tudes in the lattice. We have also shown that the occupa- 
tion of these out-of-equilibrium condensates scales with 
the number of particles in the initial Mott insulator, and 



with a proportionality constant that increases with 77. We 
followed a mean-field approach that was gauged against 
exact diagonalization results of small 2D systems, where 
qualitatively similar results were obtained within both 
approaches. 

Our results suggest that the expansion of a Mott insu- 
lator can be used to generate strongly interacting atom 
lasers in optical lattices, where the velocity of the con- 
densates (and their wavelengths) can be fully controlled 
by the parameters of the optical lattice involved in the 
setup. This study thus presents evidence for a fundamen- 
tally new phenomenon, namely, out-of-equilibrium con- 
densation in an isolated, and hence, energy-conserving 
system. As such, this study may have important impli- 
cations on current experiments with ultracold Bose gas 
in optical lattices, as well as on our understanding of 
quantum systems out of equilibrium. We note here that 
experimental setups capable of verif ying the results pre- 
sented here are currently available [l6|. As mentioned 
before, we have checked that this phenomenon also oc- 
curs in the soft-core boson case for large onsite repulsive 
interactions, in which case the dispersion relation, and 
hence the momenta at which the bosons condense, be- 
comes dependent on the value of the onsite repulsion. 
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MEAN-FIELD VS EXACT RESULTS IN ID 

In one dimension, the existence of an exact solution to 
our model of interest (see Ref. [5] in the Letter) enables 
us to compare mean-field results against exact ones in 
systems with lattices which are comparable in size with 
those realized experimentally. The ID comparison should 
however be treated with caution because the mean-field 
approximation is expected to yield meaningful results 
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only in higher dimensions. Nonetheless, in Fig.|4]we show 
several snapshots of the density profile of a ID system, 
taken at five different times. These show the mean- field 
local densities (triangles), compared against matching ex- 
act ones (circles). In both cases, 20 bosons are released 
from a trap with V = 0.01^ on a 200-site lattice. 

It is clear from the figure that despite the visible quan- 
titative differences between the mean-field evolution and 
the exact one, the two systems exhibit a qualitatively 
similar expansion. In both cases, two 'lumps' of bosons 
emerge during the melting of the Mott insulator, and 
subsequently move away from each other at the same 
constant velocities. Remarkably, the fronts of the exact 
and mean-field condensates propagate together, as can 
be verified by examining the points where the density 
goes to zero in both solutions. 
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FIG. 5: The error edens [Eq. Q] as a function of time for the 
mean-field evolution displayed in Fig. [H 



Figure [5] shows the normalized difference between the 
mean-field local densities and the exact ones, given by: 
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FIG. 4: Evolution of 20 HCBs on a 200-site ID lattice. 
The HCBs are initially trapped inside a harmonic potential 
and are then released. The figure shows the density profiles 
of the system at five different times (from top to bottom: 
r = 0, 10, 20, 30 and 40) as it was computed by exact means 
(circles) and by the mean- field approximation (triangles). 



The errors clearly increase with time. However, they do 
not increase as fast as may have naively been expected 
for a mean-field solution applied to a ID system. This is 
encouraging in terms of the relevance of the mean-field 
approach for this model when applied to higher dimen- 
sional systems. 
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MEAN-FIELD VS EXACT RESULTS IN 2D 

Here, we discuss in further detail the exact setup 
through which the analysis of the momentum-error pre- 
sented in the Letter [see Eq. (4) there] was performed. 
We also extend the analysis by considering the behavior 
of the error for two additional values of the parameter 
77 and discuss a similar analysis performed on the errors 
of the local densities [Eq. (|7J)]. This is done in order 
to strengthen the conclusions drawn in the text in that 
context. 

Since in 2D, exact diagonalization methods allow in- 
sight into the dynamics of systems with only a small num- 
ber of particles hopping on small-size lattices, the setup 
for the error analysis presented here is slightly different 
than the basic scenario discussed in the Letter where 
bosons were released from a trap whose origin was the 
center of a periodic lattice. Here, we have used a har- 
monic potential along the ^-direction centered around 
one edge of the lattice, supplemented by open boundary 
conditions in the x direction. This setup takes advan- 
tage of the parity symmetry (x — )• —x) of the expansion 
observed in the conventional setup. This method has 
been used in the time-dependent DMRG calculations in 
Ref. [8] in the Letter. In order to make the behavior 
of the different systems comparable, the characteristic 
density (see Ref. [5] in the Letter) of the different initial 
states was kept fixed for all the systems examined. This 
was done by choosing the curvature of the harmonic po- 
tential to be V = Vq/W% (in the figures, Vo = 2, and 
W x = 1, 2, or 3). In Fig. 1 in the Letter, and in Figs. [6] 
and [71 here, the behavior of the error is shown up to r = 1, 
which is the typical time it takes for the particles to reach 
the opposite edge of the lattice and bounce back. 

In the Letter, we have shown that, in general, as the 
value of 77 increases, the discrepancies between the mean- 
field and exact momentum profiles decrease. There, this 
was shown for two values of 77, namely, 77 = 0.2 and 
77 = 0.8. To complete this picture and to provide further 
corroboration of our conclusions, in Fig.[6]we display the 
behavior of the momentum-error as a function of time for 
the same small 2D lattices but with 77 = 0.4 and 77 = 0.6. 

In accord with the conclusions reported in the text, it 
can be seen in Fig. [6] that larger values of 77 yield smaller 
errors [Fig. [Bfa) vs[6fb)], but also that adding particles 
and/or sites in the transverse direction both reduce the 
errors as well. The two figures fit with the trend set by 
the 77 = 0.2 and 77 = 0.8 panels of Fig. 1 in the Letter, 
and therefore strengthen our previous argument that for 
many-particle 3D systems the mean-field approach may 
be a good approximation for the study of the dynamics 
of 3D bosons. 

Further support to our conclusions is provided by the 
behavior of the errors of the local densities. These are 
computed via Eq. (0, and depicted in Fig. [3 for the four 
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FIG. 6: (Color online) Momentum- error e m0 m [Eq. (4) in the 
Letter] as a function of time for different small 2D systems 
of size L x x L± with Nb bosons. The hopping ratios are: (a) 
77 = 0.4 and (b) 77 = 0.6. As the figures indicate, the errors 
in the 77 = 0.6 case are smaller, since larger values correspond 
to systems which are more two-dimensional. Also, as the 
number of particles in the system or the number of lattice 
sites in the transverse direction increase, the errors become 
smaller. These results are in accord with those shown in Fig. 
1 in the Letter. 



different values of 77. Here too we have found that the 
errors decrease with increasing 77, and that increasing the 
lattice size in the transverse direction and the addition 
of particles both reduce the difference between the mean- 
field and exact results. 
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FIG. 7: (Color online) The error edens [Eq. J7J)] as a function 
of time r for different small 2D systems of size L x x L± with 
Nb bosons. The hopping ratios are: (a) 77 = 0.2, (b) 77 = 0.4, 
(c) 77 = 0.6, and (d) 77 = 0.8. 
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CROSS SECTIONS 

In Fig. 2 of the Letter, we presented intensity plots of 
the density and momentum profiles of the bosonic gas 
during the expansion of a Mott insulator in 3D. In Fig. 
[H we provide the cross sections of the same density and 
momentum profiles, with cuts along the x and k x axes, 
i.e., y = z = and k y = k z = 0, respectively. 

As evident from Fig. EJa), when the initially-trapped 
bosonic gas is released, two solitonic lumps begin to form 
during the melting of the Mott insulator. Once the latter 
has disappeared, the lumps start moving away without 
a considerable change in their shape while they slowly 
broaden. The cuts across the momentum distribution 
function [Fig. E{b)] show very sharp peaks emerging at 
finite momenta. Those peaks signal condensation, as dis- 
cussed in the text, where the lowest natural orbital occu- 
pation was found to scale proportionally with the number 
of particles in the system. 
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FIG. 8: (Color online) Evolution of (a) density profiles (b) mo- 
mentum profiles (c) absolute value of the LNO wave function 
during the free expansion of a Mott insulator with W x = 40 
in a 3D lattice with L x = 700 and r\ = 0.8. The times shown 
here are r = 0, r = 84, r = 151 and r = 187. (d) Total 
number of condensed bosons (LNO occupation) in the sys- 
tem divided by the total number of bosons, as a function of 
time. 

In Fig. [8fc), we show the cross section along the x 
axis of the lowest natural orbital. The figure illustrates 
that the condensates are mainly localized in, and have 
similar shapes as, the two lobes visible in the density 
profiles. This explains why the (constant) group velocity 
of the condensates, determined by their momenta and the 
dispersion relation in the lattice (see related discussion in 



the Letter), is equal to the velocity of motion of the lobes 
appearing in the density profiles. 

The condensate fraction, i.e., the ratio between the oc- 
cupation of the LNO and the total number of bosons, is 
shown in Fig. 5(d) as a function of time. There, one 
can see that it quickly increases in the initial stages of 
the expansion, during the melting of the Mott insulator, 
and then, after the insulator has melted, it exhibits a 
slow increase during the remainder of the evolution. The 
number of condensed bosons at the kink connecting the 
two regimes is the one reported in Fig. 3(b) of the Letter. 
Interestingly, while the expansion of the Mott insulator 
in our selected geometry can only cause the emergence of 
two independent condensates traveling in opposite direc- 
tions (during their formation, the insulating slab between 
them precludes any phase correlations), the diagonaliza- 
tion of the one-particle density matrix computed within 
the mean- field approximation yields only one large eigen- 
value, i.e., mean- field merges the two LNOs into one. 
Hence, the mean-field LNO occupation reflects the total 
number of condensed bosons (that is, the sum of the oc- 
cupations of the two moving condensates). This can be 
easily verified by examining the evolution of only one half 
of the system described in the previous section. 



EFFECTS OF OPEN BOUNDARY CONDITIONS 

In order to verify that open boundary conditions in 
the transverse direction do not alter the general features 
of the systems studied in the Letter, we have tested the 
effects of open boundaries on the condensate occupation 
in the lattice as a function of the transverse linear sys- 
tem size L±. This was done in order to illustrate that 
in large systems edge effects do not destroy the general 
phenomenon described earlier. The results of this inves- 
tigation are summarized in Fig. [9] The data points rep- 
resent the condensate fraction at a fixed time (r = 18) 
after the trap has already been removed and the bosons 
were allowed to expand and condense. 

We have tested the effects of the open boundaries on 
the condensation in several 2D systems all with 100 sites 
in the x direction but with varying L±_ (here, 77 = 0.2). 
The dashed horizontal line marks the same quantity in 
the periodic (L± — >> 00 limit) case discussed in the Let- 
ter. As expected, the larger L± is, the less pronounced 
edge effects become and we can therefore conclude that 
provided L_\_ is large enough, systems with open bound- 
ary conditions (a box trap) in the transverse direction 
will exhibit the exact same features of parallel periodic 
systems. This will remain true in the time scales relevant 
for the expansion of the bosons or the formation of the 
solitonic condensates. The inset in the figure is a log-log 
scale of the data points, illustrating that the condensate 
occupation approaches the periodic value algebraically, 
where the exponent here is « —0.675. 
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FIG. 9: Effects of open boundary conditions. Condensate 
fraction (given in %) at r = 18 as a function of lattice size 
in the transverse direction, for a 2D system with L x = 100 
and rj = 0.2. As the size of the system in the transverse 
directions increases, the condensate fraction approaches the 
infinite-size (periodic) value marked by the dashed horizontal 
line. The inset shows a linear fit on a log-log scale illustrating 
the power-law behavior of the condensate fraction as it reaches 
the periodic limit. The exponent here is ~ —0.675. 



